
conf = omni_load_conf()
s = omni_read_cat(conf.cat, ncat)
IF ~exist(pvec) THEN restore,'./local/'+conf.survey+'_pvec.sav',/ver


myps,'./grs_paper/catalog_comp.eps',xsize=6.5,ysize=2.2,thick=1.5
!x.omargin = [-5.5,-5]
!y.omargin = [-0.1,-1.2]
multiplot,[3,1]
cs = 0.55

yrange = [0,1150]

plothist,s.glat,bin=0.1,/fill,fcolor='skyblue',charsize=cs,color='cyan',$
         xtit='Galactic Latitude  [deg]',xarr,yarr1,yr=yrange,$
         ytit='N per 0.1'+cgSymbol('deg')+' bin',xr=[-2,2],yst=5
cgAxis,xaxis=0,/xst,xtickformat='blank_axis'

yf1 = mpfitpeak(xarr,yarr1,A,NTERM=3)
;; cgOplot,xarr,gauss_1(xarr,A),color='green'
print,'Blue: ',[A,A[2]*2.355]


cgAxis,yaxis=0,yst=5,yr=yrange/2.,/save

ikd = where(constrain.kdist OR constrain.grsmatch OR $
            constrain.parallax OR constrain.knownd, nkd)
plothist,s[ikd].glat,bin=0.1,/over,xarr,yarr3,color='red',$
         /fill,/fline,fcolor='red',forientation=45,fspacing=0.2

yf3 = mpfitpeak(xarr,yarr3,A,NTERM=3)
print,'Red: ',[A,A[2]*2.355]

cgAxis,yaxis=0,/yst,yr=yrange/5.,charsize=cs,$
       ytit='N per bin',/save
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

iwc = where(constrain.post,nwc)
plothist,s[iwc].glat,bin=0.1,/over,xarr,yarr2

yf2 = mpfitpeak(xarr,yarr2,A,NTERM=3)
;; cgOplot,xarr,gauss_1(xarr,A),color='red'
print,'Black: ',[A,A[2]*2.355]

kstwo,s.glat,s[ikd].glat,D,prob
print,'GLAT KS: B-K: ',D,prob,100.d*(1.d - prob)
kstwo,s.glat,s[iwc].glat,D,prob
print,'GLAT KS: B-W: ',D,prob,100.d*(1.d - prob)
kstwo,s[iwc].glat,s[ikd].glat,D,prob
print,'GLAT KS: K-W: ',D,prob,100.d*(1.d - prob)



print,'RATIO2: ',max(double(yarr1))/max(yarr2)
print,'RATIO3: ',max(double(yarr1))/max(yarr3)

div = cgSymbol('div')

al_legend,/top,/right,box=0,charsize=0.9*cs,linsize=1.5,$
          ['   ','Bolocat '+div+'5',$
           'Kin. Samp. '+div+'2','W-C'],$
          color=['White','Cyan','Red','Black'],linestyle=0,thick=5/2.
al_legend,/top,/right,box=0,charsize=0.9*cs,['Bin = 0.1'+cgSymbol('deg')]

multiplot

;;====================================


plothist,alog10(s.flux),bin=0.1,/fill,fcolor='skyblue',charsize=cs,xst=5,$
         color='cyan',xarr,yarr1,yst=5,yr=yrange
;; vline,median(alog10(s.flux)),color='cyan',linestyle=3
print,'CYAN:  ',10.^(median(alog10(s.flux)))

cgAxis,yaxis=0,yst=5,yr=yrange/2.,/save

ikd = where(constrain.kdist OR constrain.grsmatch OR $
            constrain.parallax OR constrain.knownd, nkd)
plothist,alog10(s[ikd].flux),bin=0.1,/over,xarr,yarr3,color='red',$
         /fill,/fline,fcolor='red',forientation=45,fspacing=0.2
;; vline,median(alog10(s[ikd].flux)),color='red',linestyle=3
print,'RED:   ',10.^(median(alog10(s[ikd].flux)))

cgAxis,yaxis=0,/yst,yr=yrange/5.,charsize=cs,$
       ytickformat='blank_axis',/save
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

iwc = where(constrain.post)
plothist,alog10(s[iwc].flux),bin=0.1,/over,xarr,yarr2
;; vline,median(alog10(s[iwc].flux)),color='black',linestyle=3
print,'BLACK: ',10.^(median(alog10(s[iwc].flux)))

kstwo,s.flux,s[ikd].flux,D,prob
print,'FLUX KS: B-K: ',D,prob,100.d*(1.d - prob)
kstwo,s.flux,s[iwc].flux,D,prob
print,'FLUX KS: B-W: ',D,prob,100.d*(1.d - prob)
kstwo,s[iwc].flux,s[ikd].flux,D,prob
print,'FLUX KS: K-W: ',D,prob,100.d*(1.d - prob)


print,'RATIO2: ',max(double(yarr1))/max(yarr2)
print,'RATIO3: ',max(double(yarr1))/max(yarr3)

cgAxis,xaxis=0,/xlog,/save,/xst,xr=10.^!x.crange,xtickformat='exponent10',$
       xtit='Flux Density  [Jy]',charsize=cs
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'

al_legend,/top,/right,box=0,charsize=0.9*cs,linsize=1.5,$
          ['  ','Bolocat '+div+'5',$
           'Kinematic Sample '+div+'2','Well-Constrained'],$
          color=['White','Cyan','Red','Black'],linestyle=0,thick=5/2.
al_legend,/top,/right,box=0,charsize=0.9*cs,['Bin = 0.1 dex']

multiplot

;;====================================

ind = where(finite(s.rad) AND s.flux40 GT 0.)
s   = s[ind]
c   = constrain[ind]

omni_physical_constants         ; Load !PC structure with physical constants

SB = s.flux / (s.rad * s.rad) / !dpi     ; Jy / arcsec^2
SB *= (1.d-6 / 2.3504d-11)

print,m4_stat(SB)

print,'NORMALIZATION on SB: ',(1.d-6 / 2.3504d-11)/33./33./!dpi

print,m4_stat(SB)
print,m4_stat(s.flux)
print,m4_stat(s.rad)
help,SB
plothist,alog10(SB),bin=0.1,/fill,fcolor='skyblue',charsize=cs,xst=5,$
         color='cyan',xarr,yarr1,yst=5,yr=yrange,xr=[-0.15,2.65]
;; vline,median(alog10(SB)),color='cyan',linestyle=3
print,'CYAN: ',10.^(median(alog10(SB)))
print,'SB CRANGE: ',!x.crange
cgAxis,yaxis=0,yst=5,yr=yrange/2.,/save

ikd = where(c.kdist OR c.grsmatch OR $
            c.parallax OR c.knownd, nkd)
plothist,alog10(SB[ikd]),bin=0.1,/over,xarr,yarr3,color='red',$
         /fill,/fline,fcolor='red',forientation=45,fspacing=0.2
print,'RED: ',10.^(median(alog10(SB[ikd])))

cgAxis,yaxis=0,/yst,yr=yrange/5.,charsize=cs,$
       ytickformat='blank_axis',/save
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

iwc = where(c.post)
plothist,alog10(SB[iwc]),bin=0.1,/over,xarr,yarr2
;; vline,median(alog10(SB[iwc])),color='black',linestyle=3
print,'BLACK: ',10.^(median(alog10(SB[iwc])))

kstwo,SB,SB[ikd],D,prob
print,'SB KS: B-K: ',D,prob,100.d*(1.d - prob)
kstwo,SB,SB[iwc],D,prob
print,'SB KS: B-W: ',D,prob,100.d*(1.d - prob)
kstwo,SB[iwc],SB[ikd],D,prob
print,'SB KS: K-W: ',D,prob,100.d*(1.d - prob)

print,'RATIO2: ',max(double(yarr1))/max(yarr2)
print,'RATIO3: ',max(double(yarr1))/max(yarr3)


cgAxis,xaxis=0,/xlog,/save,/xst,xr=10.^!x.crange,$;xtickformat='exponent10',$
       xtit='Mean Surface Brightness  [MJy sr!u-1!n]',charsize=cs
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'

print,fix(total(finite(SB)))
print,fix(total(finite(s.flux)))
print,fix(total(finite(s.rad)))

al_legend,/top,/right,box=0,charsize=0.9*cs,linsize=1.5,$
          ['   ','Bolocat '+div+'5',$
           'Kinematic Sample '+div+'2','Well-Constrained'],$
          color=['White','Cyan','Red','Black'],linestyle=0,thick=5/2.
al_legend,/top,/right,box=0,charsize=0.9*cs,['Bin = 0.1 dex']


myps,/done,/mp

;; cgPlot,alog10(s.flux),s.rad,psym=16,symsize=0.7
;; linfit_plot,alog10(s.flux),s.rad,color='red',/legend,/top,/right
;; vline,/h,33,color='green'



print,n_elements(c)
iks = where(c.kdist OR c.grsmatch OR c.knownd OR c.parallax, nks)
print,nks
iwc = where(c.post, nwc)
print,nwc

rad2 = s.rad*s.rad

print,'Pearson:'
print,correlate(alog10(s.flux),rad2)
print,correlate(alog10(s[iks].flux),rad2[iks])
print,correlate(alog10(s[iwc].flux),rad2[iwc])

print,'Spearman:'
print,r_correlate(alog10(s.flux),rad2)
print,r_correlate(alog10(s[iks].flux),rad2[iks])
print,r_correlate(alog10(s[iwc].flux),rad2[iwc])

END
